Draft version February 1, 2008 

Preprint typeset using IAT^X style emulateapj v. 04/03/99 



R-MODES IN NEUTRON STARS WITH CRUSTS: TURBULENT SATURATION, 
SPIN-DOWN, AND CRUST MELTING 

Yanqin Wu, Christopher D. Matzner, and Phil Arras 

Canadian Institute for Theoretical Astrophysics, University of Toronto 
Draft version February 1, 2008 



o 
o 
o 

(N 



> 
m 

O 
O 
O 

^: 

Or 

6 ■ 

CO 



X 



ABSTRACT 

Rossby waves (r-modes) have been suggested as a means to regulate the spin periods of 
young or accreting neutron stars, and also to produce observable gravitational wave radiation. 
R-modes involve primarily transverse, incompressive motions of the star's fluid core. However, 
neutron stars gain crusts early in their lives: therefore, r-modes also imply shear in the fluid 
beneath the crust. We examine the criterion for this shear layer to become turbulent, and derive 
the rate of dissipation in the turbulent regime. Unlike dissipation from a viscous boundary layer, 
turbulent energy loss is nonlinear in mode energy and can therefore cause the mode to saturate 
at amplitudes typically much less than unity. This energy loss also reappears as heat below the 
crust. We study the possibility of crust melting as well as its potential implications for the spin 
evolution of low-mass X-ray binaries. Lastly, we identify some universal features of the spin 
evolution that may have observational consequences. 



1. INTRODUCTION 

The possibility of gravitational radiation from 
rapidly-rotating neutron stars has recently come un- 
der intense scrutiny. This interest has been fueled 
partly by the need to explain both the limited range 
of rotation periods in low-mass X-ray binary systems 
(LMXBs) with various accretion rates (Bildsten 1998; 
Andersson 1998) and the observed upper limit in spin 
rate of young neutron stars, and partly by the pos- 
sibility that the gravitational radiation from LMXBs 
or newborn neutron stars may be detectable by grav- 
ity observatories such as LIGO (Owen et al. 1998; 
Bildsten 1998; Andersson, Kokkotas, & Stergioulas 
1999). For neutron stars with crusts, there are two 
potential sources of gravitational radiation: a mass 
quadrupole originating in the solid crust of the star 
(Bildsten 1998), or a current quadrupole from Rossby 
waves, which become overstable under the influence 
of their own gravitational radiation (Andersson 1998; 
Friedman & Morsink 1998). 

The outer layer of a neutron star solidifies if the 
typical energy of Coloumb interactions between nu- 
clei, Z 2 e 2 /r, exceeds the thermal energy, kf,T, by 
a critical factor 172 (Farouki & Hamaguchi 1993). 
Here, Z and r are the mean charge and the mean 
spacing of nuclei, respectively. The crystallization 
(also the melting) temperature is therefore 
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where we have scaled Z and A (mean atomic weight) 
with values suitable for the bottom of the LMXB 



crust (Haensel & Zdunik 1990). Young neutron stars 
(Negele & Vautherin 1973), and LMXBs whose crust 
have been sufficiently heated, may have higher val- 
ues of Z and therefore higher T me it Nucleation of a 
crystalline crust begins at around T me i t , and crystal- 
lization accelerates exponentially as the temperature 
is decreased further (de Blasio 1995). All but the 
hottest neutron stars are expected to possess a solid 
crust of thickness ~ 1 km. Our investigation contin- 
ues along the lines of Bildsten & Ushomirsky (2000; 
hereafter, BU) to explore the implications of such a 
solid crust for r-modes. We also consider the fate of 
the crust during r-mode instability. 

We shall focus on the fastest growing r-mode with 
spherical indexes i = m = 2, whose rotating frame 
frequency is u = 2S7 s /3. Here £l s is the spin frequency 
of the neutron star; we denote its angular counterpart 
Q s /2tt. In the fluid core, the r-mode is largely 
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horizontal, and the modulus of its displacement vec- 
tor is (Owen et al. 1998) 
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where 9 is the co-latitude, r the distance from the cen- 
ter and R the stellar radius. The amplitude of the ve- 
locity perturbation is related to the displacement by 
\v\ = o>|£|. If the neutron star crust is perfectly rigid, 
the r-mode produces a periodic rubbing at the fluid- 
solid boundary with velocity \Av\ ~ How- 
ever, Levin & Ushomirsky (2000) have argued that 
the crust of a neutron star is not perfectly rigid, and 
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across the core-crust boundary. The r-mode periodi- 
cally rubs the fluid core against the solid crust with 
this velocity. A viscous boundary layer develops for 
small velocity jump, and a turbulent boundary layer 
sets in when the jump increases above some critical 
value. 

The effect of a viscous boundary layer on r-modes 
was first analyzed by BU. We reiterate some of their 
results here. Ignoring the Coriolis force, the thickness 
of this layer is given by the diffusion length during one 
mode period (Stokes 1851), 



(3) 



where v is the molecular viscosity at the base of the 
crust, v ns 1.8xl0 4 T 8 " 2 cm s for the chosen density 
(Flowers & Itoh 1979; Cutler & Lindblom 1987), and 
Tg stands for T/10 8 K. The rate of energy dissipation 
per unit area is equal to the mean relative kinetic 
energy contained within this layer (half of p\Av\ 2 /2 
times 5) multiplied by the mode frequency: 
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it will participate in the lateral motion of the mode. 
To account for this, we write |Auj = rj\v\, with r] = 1 
in the limit where the shear force in the crust much 
exceeds the local Coriolis force. 

BU considered energy dissipation in a thin vis- 
cous boundary layer at the above core-crust interface. 
Dissipation due to "molecular" viscosity is more in- 
tense in this laminar boundary layer than in the bulk 
flow caused by the mode. Therefore, the presence of 
the boundary layer increases the stellar rotation rate 
above which the mode becomes unstable, relative to 
the case when the crust is absent (Owen et al. 1998). 
The mode is pumped by gravitational radiation and 
damped by viscosity at rates that both scale as a 2 . 
So a laminar boundary layer cannot provide a satu- 
ration mechanism for an unstable r-mode. 

In contrast, a turbulent boundary layer provides 
nonlinear dissipation that leads to r-mode saturation. 
A turbulent boundary layer occurs once the mode 
amplitude grows above a critical value. It removes 
kinetic energy from the mode with a rate that is cu- 
bic in a. As a grows, turbulent dissipation rapidly 
catches up with the amplification of the mode en- 
ergy due to gravitational back-reaction. We find this 
mechanism typically saturates the mode amplitude a 
at values much less than unity, i.e., well below what 
has been assumed in the literature (Levin 1999; Owen 
et al. 1998) 

The kinetic energy of the r-mode is converted into 
heat in the thin turbulent boundary layer. 1 If con- 
duction and local neutrino emission cannot carry the 
heat away sufficiently quickly, the local temperature 
may increase above the melting temperature (eq. [0) 
and the crust will begin to melt. In §||, we calculate 
the thermal profile induced by boundary-layer heat- 
ing at the core-crust interface. 

Finally, in §||, we consider the implications of tur- 
bulent saturation and local heating on the thermal 
and spin evolution of LMXBs as well as new-born 
neutron stars. We conclude in §[|. 

Throughout this paper we employ fiducial 
model an n = 1 polytrope of mass M = 1.4M and 
radius R = 12.53 km (as in Lindblom, Mendell, & 
Owen 1999). The bottom of the crust is assumed to 
have density p = 1.5 x 10 14 gem" 3 (Brown 2000). 

2. VISCOUS AND TURBULENT BOUNDARY LAYERS 

Going outward from the fluid core to the solid crust, 
the shear modulus of the material changes from zero 
to a large (but finite) value. This implies a discon- 
tinuity Av in the horizontal velocity of the r-mode 

x The heat input due to the viscous boundary layer is insignificant except when the r-mode is marginally unstable. 

2 This expression for E assumes that the r-mode satisfies equation (^) in the whole star; strictly speaking, an impenetrable crust 
reduces the mode energy by a factor of ~ 2. We ignore this correction throughout the paper. 

3 Note that if we define a Reynolds number using the width of the turbulent layer (as is done in a pipe flow), as opposed to using 
A£, we find K cr it ~ 500 at the transition to turbulence (Landau & Lifshitz 1959). 
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This yields a global energy dissipation rate of 

E vU « -0.03 rVT^, (5) 



where fkHz = v 8 jl kHz, and E refers to the energy in 
the mode (Owen et al. 1998), 



E = -a 2 n 2 MR 2 J, 
2 s 
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with J = 0.016. 2 BU found that the VISCOUS 
boundary layer stabilizes r-modes against gravita- 
tional back-reaction excitation at low spin frequencies 
and low stellar temperatures. 

When an unstable r-mode acquires sufficiently large 
amplitude, the laminar boundary layer will make 
a transition to turbulence. This requires the local 
Reynolds number to exceed a critical value 
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where the amplitude of the relative displacement 
|A£| = | A?; I/a; and the critical Reynolds number 
is determined experimentally to be 5R cr it ~ 2 x 10 5 
(Jensen, Sumer, & Fredsoe 1989). 3 
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Combining the Reynolds number criterion (|7|) with 
the definition of |A£|, we find that the onset of tur- 
bulence at the equator requires a mode amplitude a 
greater than a critical value a cr it: 
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which corresponds to a thickness of the boundary 
layer d = 31cmT g ~ z^,, (see Appendix |A|) . At other 
latitudes, turbulence sets in at a larger mode ampli- 
tude. 

Turbulence implies that eddies, rather than viscos- 
ity, advect momentum toward the crust. The crust 
feels a drag force very nearly in phase with Av with 
a magnitude (eq. |A6[|) 



Drag = C D -p\Av\ 2 /2, 



(9) 



where the drag coefficient Cd is measured to be 
Cd ~ 5x 10 -3 (Jensen et al. 1989) at the onset of tur- 
bulence and varies logarithmically with the Reynolds 
number beyond the critical point, as is discussed in 
Appendix [A]. Kinetic energy is removed from the r- 
mode with a time-averaged rate per unit area of 
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where time-averaging gives the factor 1/2 as the drag 
is in phase with Av. Integrating over the sphere as- 
suming a constant Cd, we find the total dissipation 
rate 

37F ^2^_| A „. |3 
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where \Avi\ = (5/167r) 1 / 2 a£l s R is the velocity jump 
at the equator. In general, we adopt for Cd the value 
at the equator, where the Reynolds number is the 
largest. 

Turbulent dissipation increases with a more rapidly 
than does the pumping of the mode by gravitational 
back-reaction, and this leads to saturation. Recall 
that the mode gains energy due to gravitational wave 
radiation at the rate (cf. Owen et al. 1998) 
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Balancing this energy gain with the nonlinear energy 
drain from turbulence, we find a saturation amplitude 



a sat ~ 3.5 x 10 3 i] 
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Note that a sa t is independent of temperature up to 
the logarithmic dependence of Cd on K. Typically 
«sat > «crit for an unstable r-mode, in which case 
the r-mode grows rapidly until it reaches a steady 



state with a = a sat . The scaling a sat oc i]~ 3u kUz 
seems to imply that the saturation amplitude will in- 
crease enormously (at a given spin frequency) if the 
crust's rigidity parameter rj is reduced from unity to 
0.1. Newborn, rapidly spinning neutron stars will 
therefore spin down at a rate that is very sensitive to rj 
(and may be saturated by other means if eq. [O pre- 
dicts a sa t > 1). In contrast, the spins of LMXBs are 
confined within a narrow frequency range that varies 
with T] (as we will show in §|j). The spin frequency at 
which they first become unstable scales as v s oc 77 4 / 11 
(equating eqs.j| with |i~2|| ). Taking this change of 
spin frequency into account, we find that the max- 
imum r-mode amplitude in LMXBs varies as r]~ 1 ' 2 . 
Low saturation amplitudes are expected in LMXBs 
even for rj ~ 0.1. 

We return now to comment on a few issues relevant 
to turbulence onset. The first is related to the local 
stratification. We find the Richardson number 
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at a = ct sa t for a Brunt- Vaisaia frequency of N ~ 
500 rad s -1 (Reisenegger & Goldreich 1992). This in- 
dicates that the stable stratification at the core-crust 
boundary does not prevent the onset of turbulence. 

Notice that the Reynolds number (eq. ^) can also 
be written as 5ft = 4(|A£|/J) 2 ~ 2C D (d/5) 2 , where 
the thickness of the turbulent layer d is related to 
|A£| as d ~ (C D /2) 1 / 2 |A^| (see Appendix |Sj). The 
condition > 3? cr i t implies d S> 5. Viscosity is unim- 
portant in most of the turbulent region. Moreover, 
the turn-over time of the energy-bearing (also the 
largest) eddies in the turbulent region is much shorter 
than the mode period, \Av\/d ^> u. This allows for 
a well-developed turbulent cascade. 

The interaction between shear turbulence and stel- 
lar rotation, the possibility of equipartition magnetic 
fields in the turbulent layer, and the possible super- 
fluid nature of the material, are potentially important 
questions that lie beyond the scope of this paper. 

3. TEMPERATURE PROFILE FROM BOUNDARY 
LAYER HEATING 

Turbulence converts the kinetic energy of the r- 
mode into heat and deposits it in a thin region of 
width d ~ 100 cmr/- 2 z^ Hz below the crust. If thermal 
conduction and local neutrino emission cannot carry 
away the heat efficiently, the temperature in this layer 
may be driven up to past the melting temperature of 
the crust (eq. ||). As our saturation mechanism (§||) 
depends on the presence of such a solid crust, we wish 
to establish the temperature profile of the boundary 
layer. We assume a steady state solution because the 
conduction time over the regions of interest to us is 
much shorter than the evolution time. 
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The conduction of heat away from the bound- 
ary layer requires a temperature profile that declines 
in each direction. This, in turn, implies that the 
rate of local neutrino cooling decreases with distance 
away from the boundary layer. We assume that the 
scale length of the temperature perturbation is much 
shorter than the density scale height, and therefore 
much smaller than the stellar radius as well. This as- 
sumption, whose validity must be checked a posteri- 
ori, allows us to consider only planar heat conduction 
in a region of constant density. 

The following toy model gives insight into the full 
solution. Let us consider a piecewise linear tempera- 
ture profile, one that reaches a temperature T b i within 
the boundary layer (of width d), and falls to zero over 
lengths l core and l crus t to the inside and outside, re- 
spectively. The two lengths differ because the ther- 
mal conductivity takes different values in the two di- 
rections. 

In a state of equilibrium, any energy not radiated as 
neutrinos within the boundary layer must be carried 
away by conductive fluxes: 
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where neutrino emissivity due to modified Urea reac- 
tions is e v ~ £qT s erg cm" 3 s -1 with eo = 7.4 x 10 -52 
(Friman & Maxwell 1979). The heat flux leaving the 
boundary layer in each direction is 

and -Fcrust = K crust 7^-, (16) 



lc. 



respectively. The thermal conductivity of the crys- 
talline crust is dominated by electron-phonon scatter- 
ing and scales with the local temperature inversely, 

Kcrust = ^crustTg" 1 RS lO 20 ^" 1 erg Cm" 1 S" 1 K" 1 , 

while that of the liquid core is dominated by electron- 
proton and electron-neutron scattering with the 
same temperature dependence, K core = RcoreT^ 1 m 
lCFTg^ergcm-is^K" 1 (Flowers & Itoh 1979). 

Another consequence of thermal equilibrium is that 
the heat conducted into the crust or core must be ra- 
diated there in neutrinos. With our approximation of 
a linear temperature profile, this implies F ~ 
in each direction with the conductive lengths I de- 
fined in equation (|l~6|). A more accurate calculation 
in which no temperature profile is assumed a priori 
(Appendix [Bj) gives F = eqIT^/A. Combining this 
with equation (|16|), we find the conductive lengths to 
be 
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2.4 x 10 5 cm (l0 9 K/T b i 
7.4 x 10 3 cm (lO 9 K/T b i) 4 . (17) 



As long as Ty > 5 x 10 8 K, / C rust < 'core "C R and our 
assumption of a thin heated region is justified. The 
heating becomes increasingly local as the boundary 
layer rises in temperature: when more heat is de- 
posited in the boundary layer, a larger temperature 
gradient is required to conduct the heat outward lead- 
ing to a smaller conduction length. The assumption 
of an isothermal star is thus not justified in the pres- 



ence of any significant localized heating. In §4.1 and 
§ 4.2 we shall compare the implications of the isother- 
mal approximation with those of localized tempera- 
ture profile for neutron star spin evolution. 

Substituting equations ( fl6| ) and ([17]) into equation 
(|l~5|), adopting a = a sat and evaluating dEturb/dA 
at the equator, we find a quadratic equation for T bl , 
whose solution gives 

T bl = 1.5 x 10 10 K rT 1/2 *g£ 8 (v / TTb 2 - b) ^ , 

(18) 

where b = 7.4 x lO -3 ^ 4 ^^^ 2 and Cd is taken to 
be 0.005. Neutrino losses from the boundary layer 
itself dominate over the conductive fluxes for b < 1, 
or u s > 650 Hz ry 8//23 . For a given mode amplitude, 
the boundary layer temperature decreases towards 
the poles. 

We find that the flux going to the core is roughly 
(ftcore/^crust) 1 / 2 ~ 30 larger than that conducted to 
the crust. However, if the above equilibrium tempera- 
ture exceeds T me i t , the crust will begin to melt. Once 
melting begins, we expect that the energy input from 
the boundary layer will be used almost entirely to 
melt the crust. The temperature in this layer will be 
kept very close to T me it by the continual inflow of cool, 
freshly molten crust, as the boundary layer moves up- 
ward following the bottom of the receding crust. The 
time it takes to completely melt the crust can be esti- 
mated to be a few days (the ratio 10 erg/E gr where 
10 48 erg is roughly the chemical binding energy of the 
crust). 

4. IMPLICATIONS FOR NEUTRON STAR SPIN AND 
TEMPERATURE EVOLUTION 

We now consider the implications of turbulent sat- 
uration for the evolution of young neutron stars and 
LMXBs in the (£l s ,T) plane. In this section we as- 
sume the crust is present and is not melted, an as- 
sumption we examine in § |4.2| . Once the neutron star 
is pushed into the instability region, either by cool- 
ing or by accretion, its r-mode quickly grows to its 
saturated value as given by equation (|l~3|). Inserting 
it into equation (C2) we find that the time to spin 
down to a frequency f^Hz is fairly independent of the 
initial spin frequency of the star, and is given by 
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The time to spin the star down to, say, v s = 500 Hz 
is roughly lO 3 // 6 yr. As a sa t decreases sharply with 
decreasing frequency, most of the spin-down time is 
spent near the lowest frequencies. 

The above conclusion on the evolution of the spin 
frequency may need modification if the crust can be 
melted by the r-mode heating. This prompts us to 
study the temperature evolution of the star. This we 
do by first adopting the isothermal approximation for 
the whole star. It applies in the limit of infinite ther- 
mal conductivity or large-scale heating. We then ap- 
ply our results from §|| to the realistic case in which 
the conductivity is finite and r-mode heating is local. 

Evolution in the unrealistic isothermal case can be 
easily integrated using equations in Appendix |C[ It 
is useful for comparison with previous work. More 
importantly, insights gained from the isothermal case 
can be applied to the realistic case. 

4.1. Isothermal Approximation 

After the r-mode reaches its saturation amplitude, 
it provides the neutron star with a source of heat 
from turbulent dissipation. The neutron star heats 
up quickly until it reaches a state in which the r-mode 
heating is balanced by the neutrino losses. With time 
the temperature will decrease slowly because heating 
by r-mode decreases with the spin frequency. We set 
the r-mode heating term equal to the neutrino cool- 
ing term in equation ( |C3| ) to find the following scaling 
between the temperature and the spin rate of the star, 
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where we take a = a sa t and ignore the logarithmic de- 
pendence of Cd on temperature by using Cd = 0.005. 
This relation gives us the thermal equilibrium spin- 
down line. It intersects the instability curve at the 
terminal frequency vt — 520 Hz r/ 0,35 . All neutron 
stars entering the instability region with an initial 
spin frequency v s > vt will converge onto the equi- 
librium line and exit at vt- The ones that enter at 
v s < ut find their spin rates hardly affected by r- 
mode instability. This is because their spin down 
time (eq. ]l9|]) much exceeds their neutrino cooling 
time. 

Figures [l] and § exhibit the results of our evolu- 
tionary calculations with a variety of mass accretion 
rates in LMXBs and a range of initial spin rates for 
young neutron stars. The evolution equations we 
adopted for the integrations are listed in Appendix 
[C]. They are similar to those derived in Owen et al. 
(1998) apart from the additional terms that arise 
from boundary layer dissipation. We consider only 
rj = 1 here. We find that, up to logarithmic depen- 
dence of Cd on the Reynolds number, the numerical 
results confirm equations (|l~9|) and (|20|), 



Note that the presence of an equilibrium spin-down 
line is a universal feature independent of the satura- 
tion mechanism. If r-modes are saturated to constant 
values, for instance, we find the equilibrium temper- 
ature T oc v s . 

Levin (1999) was the first to point out that LMXBs 
undergo limit cycles in which the angular momentum 
accreted over the whole cycle is radiated by grav- 
itational waves when the r-mode is unstable. The 
r-mode phase lasts a time that depends on the satu- 
ration mechanism and is relatively insensitive to the 
initial spin frequency at which the r-mode is destabi- 
lized. Equation (19) yields the spin-down time at the 
terminal frequency vt to be ~ 10 3 r/ 10 / 31 yr. Com- 
pared to an accretion spin-up time of ~ 10 7 yr at 
the Eddington accretion rate, we find a maximum 
duty cycle of ~ 10~ 4 during which time the gravi- 
tational wave emitted by the r-mode may be poten- 
tially detectable, with the duty cycle decreasing at 
lower accretion rate. This duty cycle is much higher 
than that found by previous investigators (e.g., Levin 
1999) mostly due to our much reduced saturation am- 
plitude. 

4.2. Realistic Evolution Accounting for Local 
Heating 

For simplicity, we consider only the evolution of 
LMXBs in this section. Numerical integration for the 
evolution tracks in the realistic case is much more 
difficult than in the isothermal case. Fortunately, 
they are not necessary. As is evident from Figure |], 
given initial conditions, the equilibrium line and the 
instability curve completely determine the evolution 
tracks. 

The r-mode instability curve in the low tempera- 
ture range depends only on the boundary layer tem- 
perature because damping from the viscous bound- 
ary layer dominates over other damping mechanisms. 
The equilibrium spin-down line is also given in terms 
of the boundary layer temperature (eq. |jl8[ ) . So we 
can simplify the problem in hand by considering only 
the temperature at the boundary layer. 

There are two main differences between the realistic 
and the isothermal lower terminal frequency 

and an enhanced heating in the boundary layer. 

The new equilibrium spin-down line (eq. fl8]l ) in- 
tersects the instability curve at a terminal frequency 



4907? - 35 Hz 
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that is lower than in the isothermal case. It takes 
slightly longer to spin the star down to this lower 
value of vt- 

The r-mode heats the boundary layer to a higher 
temperature than under the isothermal approxima- 
tion. This results from the fact that, unlike in the 
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Fig. 1. — Evolutionary tracks for LMXBs and young neutron stars under the isothermal approximation, with arrows indicating 
the direction. The thick lines extending up from the bottom enclose the r-mode instability region. The dashed loops labeled 8 — 11 
are LMXB tracks for accretion rates M/M = 10 -8 — 10 _11 yr _1 . The long- dashed- short- dashed curves are young neutron star 
tracks with a range of initial frequencies. Note that all tracks that enter the instability region above vt — 520 Hz spin down along 
the equilibrium line (eq. po|) and exit at v s — vt (labeled as c). Below this frequency the cooling time is much shorter than the 
spin-down time (eq. Symbols a and b mark the place where two specific tracks first hit the equilibrium line. See Figure |^ 

for more details on these two tracks. 





io- 6 io- 5 io- 4 io- 3 io- 2 0.1 1 10 1 10 z 10 3 



t(yr) 

Fig. 2. — R-mode amplitude a, stellar spin frequency u s and temperature T are plotted here as functions of time along two evo- 
lutionary tracks: a young neutron star with an initial frequency u s — 1kHz {dashed lines) and an LMXB with M/M = 10 -8 yr _1 
[solid lines). The r-mode growth rate in the case of the LMXB does not become appreciable until accretion has pushed the star 
sufficiently far into the instability region. Both tracks exit the instability region at the same spin frequency (vt) after spending 
roughly the same amount of time (eq. Jis|). Symbols a and b mark the points when r-mode heating is first balanced by neutrino 
cooling and the stars begin to evolve along the equilibrium spin down line. This occurs after r-modes have reached the saturation 
amplitudes. 
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Fig. 3. — Limit cycles for the same four LMXBs in Figure |l| are plotted as dotted loops against the temperature at the boundary 
layer. The heavy dashed curves delineate the instability regions and the slanted solid lines show the equilibrium spin-down curves 
(eq. for two different values of r\. The large circles denote vt for each case. The light dashed lines are drawn for comparison 

- they are equilibrium spin-down curves under the isothermal approximation. We also show the melting line (dot-dashed) with 
T melt = 6 x 10 9 K. 



isothermal case, there is now only a shell of thickness 
(^core + ^cmst)/4 + d which receives heat and cools by 
neutrinos (^3|). 

We present the LMXB limit cycles for the two cases 
7] = 1 and 77 = 0.1 in Figure ||. One observational con- 
sequence suggests itself in the figure: all LMXBs are 
expected to have spin frequencies falling in the narrow 
range of 490 - 700 Hz when 77 = 1, and 220 - 300 Hz 
when 77 = 0.1, independent of their accretion his- 
tory. This continues the idea proposed by Bildsten 
(1998) and Andersson et al. (1999) that r-modes may 
be instrumental in halting the LMXB spin-up. The 
resemblance between the above frequency range for 
77 = 0.1 and the observed LMXB spin rates (van der 
Klis 2000) is intriguing. However, it may be difficult 
to explain the fastest millisecond pulsars (period of 
1.5 ms) using the same value of 77. 

The above results may need modification if the 
crust melts during the r-mode evolution. Melting the 
crust is more likely than is suggested by the isother- 
mal approximation. As Figure |3] demonstrates, the 
crust is most likely to melt if 77 is large and if the 
accretion rate is low. We define a melting frequency, 
I'm, at which the equilibrium line intersects the melt- 
ing line T = T melt . We find for T melt = 6 x 10 9 K, 
v m = 660 Hz for 77 = 1 and 315 Hz for 77 = 0.1. We 
speculate below on some possible evolutionary conse- 
quences when crust melting is taken into account. 



4.3. Possible Outcomes Including Crust Melting and 

Forming 

The rotational energy (~ 10 51 f^ Hz erg) far exceeds 
the chemical binding energy of the crust (~ 10 48 erg) 
for spin rates of interest: crust melting is not inhib- 
ited by lack of energy. Moreover, as discussed in §||, 
the temperature at the bottom of the crust may be 
driven above the melting temperature by turbulent 
heating; this causes melting. For our fiducial param- 
eters, this occurs for slowly accreting LMXBs if 77 = 1 
(Fig. i- 

After the crust is melted, some other nonlinear dis- 
sipation must cause the r-mode to saturate. It is 
beyond the scope of this paper to predict the evolu- 
tion when the crust is molten. However, we can be 
confident that the star will cool and spin down suffi- 
ciently to form a crust again. Once this occurs, the 
subsequent evolution is again within our jurisdiction, 
and it depends on the stellar spin frequency (f e ) at 
the moment the crust forms. We discuss five repre- 
sentative evolution scenarios that begin with different 
values of v e , labeled A through E in Figure ^j. 

In case A, the stellar temperature drops below the 
melting temperature while the star is still spinning so 
fast [y e > v m ) that the boundary layer temperature 
would exceed T me it if a crust did form. Therefore, 
the material that could form a crust is kept hot (by 
some means) with T ~ T me i t . This implies that the 
formation of an appreciable crust is delayed until the 
rotation has slowed to v s = z/ m , where v m is the fre- 
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Fig. 4. — Possible outcomes of spin evolution when we consider crust melting and forming. The horizontal axis represents the 
temperature of the boundary layer as long as this is less than the melting temperature (the dot-dashed line, we takeT mc it = 6xl0 9 K 
as is appropriate for nuclei with Z = 20); for higher temperatures this axis represents the temperature of the fluid core. The two 
dashed lines are instability curves that are valid in the presence (to the left of the melting line) or absence (to the right) of a crust, 
respectively. Equation ( |l8| ) gives the equilibrium, saturated spin-down behavior (solid line) when local heating is considered. We 
demonstrate the evolution (dotted lines with arrows) for five scenarios in which the crust forms at different initial spin frequencies. 
We find their final rotation rates reside in the range marked as the 'exit range'. In this figure, we have taken rj = 1 for illustrative 
purposes. 



quency at which the equilibrium line hits the melting 
line. After this point the star will spin down with 
a saturated r-mode (eq. [fL9f ) and cool in a state of 
thermal equilibrium (eq. |J18j|); this behavior is repre- 
sented by the solid line in Figure ^. The star leaves 
the instability region with a specific (terminal) spin 
frequency when the mode is no longer unstable. 

Case B represents the outcome if the initial spin 
frequency v e is low enough to allow crust formation, 
yet higher than the terminal spin frequency. A crust 
forms and the star evolves towards the equilibrium 
spin-down line, with the the cooling of the bound- 
ary layer regulated by that of the whole star. The 
subsequent evolution is identical to case A. 

Case C enters with v s < i>t- Its spin-down time 
is longer than the cooling time and it exits the r- 
mode instability curve without being significantly 
spun down. Case D has a similar outcome: r-mode 
instability is quenched entirely once the crust forms. 

Case E hits the instability curve to the right of the 
melting line. It evolves along the critical stability 
curve until the melting line, after which the r-mode 
becomes stable. This is the only scenario in which we 
can predict the evolution of an r-mode in the absence 
of a crust. 

From these considerations we conclude that the fi- 
nal rotation rates of young neutron stars, and of those 
LMXBs capable of melting their crusts, depend on 



their spin evolution while the crust is absent. In the 
case that r\ = 1 and T me i t = 6 x 10 9 K, the final spin 
frequencies should lie in the range of 50 to 480 Hz, 
where the upper limit is the terminal frequency vt- 
The chosen value of r] should be considered uncertain, 
but recent calculations (Levin & Ushomirsky 2000) 
point to rj ~ 0.1. Also, a minimum value of r/ ~ 0.003 
can be derived from equation (21) by assuming that 
i>t is larger than the frequency of the fastest known 
young neutron star (~ 16ms; Marshall et al. 1998). 
Moreover, young neutron stars, or LMXBs that have 
been heated above the temperature necessary for nu- 
clear burning, are more likely to be in a state of chem- 
ical equilibrium at the base of their crusts; this im- 
plies higher values of Z, thus higher melting temper- 
atures. However, the above discussion demonstrates 
that this would only raise the lower bound for the 
possible range of final rotation rates without chang- 
ing the upper bound (fy). 

5. CONCLUSIONS 

R-modes in the fluid core of a neutron star gener- 
ate velocity shear below the solid crust. The primary 
conclusion of this paper is that this velocity shear can 
drive turbulence which in turns limits the amplitude 
that an overstable r-mode can obtain. This satura- 
tion amplitude is easily calculable. It rises steeply 
with mode frequency and depends on the tempera- 
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ture at the turbulent boundary layer only logarithmi- 
cally. For most spin frequencies of interest, we find 
this amplitude falls well below unity. 

We also find that heating due to turbulent dissipa- 
tion in the boundary layer leads to a localized temper- 
ature enhancement relative to the surrounding mate- 
rial, and the heat is removed by neutrinos from both 
the boundary layer itself and from the nearby regions 
to which the heat is conducted. It is erroneous to 
assume that the whole star remains isothermal in the 
presence of a localized heat source. The boundary 
layer temperature may exceed the melting tempera- 
ture of the crust. 

When studying the spin and temperature evolu- 
tion of neutron stars undergoing r-mode instability, 
we find the following remarkable behaviors: because 
of the strong dependence of neutrino cooling on the 
temperature, the temperature of the boundary layer 
rapidly reaches an equilibrium value that depends on 
the mode amplitude and spin frequency, and hence, 
for a saturated mode, only on the stellar spin. This 
defines an equilibrium spin-down line along which all 
neutron stars evolve. We define the spin frequency 
at which this line intersects the instability line as vt, 
the terminal frequency. This frequency is a signif- 
icant feature of the spin evolution. We expect all 
neutron stars with crusts to exit the r-mode instabil- 
ity range with v s = ut, as long as their initial spin 
frequencies are higher than vt- Neutron stars born 
with v s < vt will never experience significant spin- 
down. If the crust can be melted and reformed, the 
terminal frequency defines the upper limit to the final 
spin rates of these neutron stars. We also establish a 
spin-down law for which the time required to spin the 
star down to vt is roughly independent of the initial 



spin rate. These features are expected to persist for 
other saturation mechanisms. 

This work has relied on several assumptions and 
approximations which should be noted. We have 
ignored any dissipation that may occur within the 
solid crust itself, the possibility that the crust may 
be melted by mechanical strain or heat from nuclear 
reactions, and additionally any nonlinear dissipation 
caused by coupling between modes in the bulk of the 
fluid core (Schenk et al. 2000). We have ignored the 
possibility that gravitational radiation may be dom- 
inated by a permanent mass quadrupole within the 
crust. We have not considered the effects of super- 
fluidity in the material below the crust. The detailed 
value of the boundary layer temperature depends on 
the exact value of the thermal conductivity as well 
as on the parameter rj that describes the crust mo- 
tion; these we consider relatively uncertain. More- 
over, melting of the crust is sensitive to the compo- 
sition at its base, and this may differ between young 
and accreting neutron stars. Finally, our discussions 



in §4.3 do not include the possibility that chemical 
composition in the crust may change when melting 
occurs. 
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APPENDIX 

TURBULENT BOUNDARY LAYER 

In this appendix, we give a short description of the turbulent boundary layer caused by shear flows. We 
provide scalings for the thickness of the layer, the velocity profile inside the layer and the drag force acting on 
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the wall. These scalings are confirmed by numerical experiments and we use them in the main text to derive 
the energy dissipation rate due to shear turbulence. 

Let Av be the oscillatory rubbing velocity between the neutron star's fluid core and its solid crust. The 
horizontal displacement is = Av/tu. When the boundary layer becomes turbulent, we can describe it with 
results from a similar situation: the turbulent boundary layer generated by a steady flow along a semi-infinite 
plane. Our discussion therefore follows closely that of Landau & Lifshitz (1959, §42-44). In this case, A£ 
can be understood as the distance away from the edge of the plane, and Av as the velocity of the shear flow. 
Although the steady solution is expected to reproduce the correct scalings for the boundary layer structure, it 
is time-independent and therefore cannot predict the phase shift between the shear velocity and the drag force. 
We obtain this phase shift from the experimental data of Jensen et al. (1989). 

The thickness of the turbulent boundary layer (d) grows linearly with the distance away from the edge of the 
plane (A£). We write 

d = AA£, (Al) 

where A is a small dimensionless number that describes the expansion of the turbulent layer relative to the 
horizontal displacement away from the edge. If we define v * to be the root mean squared fluctuating velocity 
in the turbulent layer, we find A = v*/Av 

Within the boundary layer, the time-averaged velocity u at distance y away from the plane obeys a logarithmic 
profile, 

„(, ) = ^ ln f£jMV (A2) 



ci V v 

where c\ ~ 0.4 and C2 — 7.69 are experimentally determined constants. At y = d, we should recover the 
flow velocity, so Av = v^c^ 1 In^v^d/u). Combining this expression with equation (Al), we find the following 
solution for A, 

j = In (c 2 A 2 K) , (A3) 
where we define an experimentally controlled Reynolds number for the flow by 

Rss^. (A4) 

A rough scaling for A can be found in the limit of very large Reynolds number where the solution to equation 
S is 

The fact that A is a slowly varying function of the Reynolds number helps simplifying the calculation for energy 
dissipation rate. 

The drag force per unit area on the surface due to the fluctuating velocity field is 

Drag = pvl = \ 2 pAv 2 = ^C D pAv 2 , (A6) 

which identifies the drag coefficient in the turbulent regime to be Cd = 2A 2 ~ 2[ci/ ln(c25R)] 2 , also a slowly 
decreasing function of the Reynolds number. Jensen et al. (1989) show that in the fully turbulent regime, the 
drag force is nearly in phase with the shear velocity. 

The scalings Cd = 23ft - 1//2 in the laminar regime (Landau & Lifshitz 1959) and Cd = 2 A 2 in the turbulent 
regime compare very well with the experimental determinations of Co by Jensen et al. (1989) over a large 
range of Reynolds number. Furthermore, if we estimate the critical Reynolds number for the onset of turbulence 
by equaling the two Cd expressions, we find Kcrit = 1.6 x 10 5 , in good agreement with the experiment. At 

= Wtcriti ^ = 0-05 and Cd = 0.005. Beyond the onset of turbulence, Cd decreases logarithmically with the 
Reynolds number. 

EQUILIBRIUM TEMPERATURE PROFILE AROUND A LOCALIZED SOURCE OF HEAT 

We assume thermal equilibrium has been reached. On either side of a localized heat source, the heat conduc- 
tion and energy conservation equations read, 

F = -kVT, 

V-F = -e v , (Bl) 
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where k is the thermal conductivity and scales inversely with T in the neutron star interior, while e„ is the 
neutrino emissivity with t v oc T 8 for modified Urea reactions. We define two constants kq = kT and eo = e u /T s . 
Under the planar approximation, we manipulate the above equations to find 

dlnT d 2 lnT = £0 T s d\nT 

dr dr 2 kq dr 



where r is the radius. Integrating both sides, we find 

'dlnT\ 2 
dr 



4At 



(B3) 



where a and b are two arbitrary points within the integration range. Here, we take point a to be at the boundary 
layer, which is the source of heat, and point b to be at many conduction lengths away from the boundary layer 
so that the flux and temperature at point b may be regarded as negligible. The flux and temperature at the 
boundary are then related by 

F=i(6 K ) 1 / 2 r b 4 1 . (B4) 
If we define a conduction length I = kT^i/F = kq/F, we find F = eo a result we adopt in §|3[ 

EVOLUTION EQUATIONS 

We have modified the phenomenological model of Owen et al. (1998) to calculate the effect of the r-mode 
on the spin and temperature of the neutron star in the isothermal approximation. The evolution equations for 
mode amplitude a, spin frequency O s , and stellar temperature T are 

d = a a l - a 2 Q 

l + a 2 Q' v ' 
• _ 2Q S a 2 Q 1 /4\V 2 M 

n °--— tt^q + tUJ M nd > (C2) 

C T T 9 f 9 = -L U>9 T$ + e a (M/M)(GM 2 /2R) + 2£/r heat , (C3) 

where Q = 0.094, / = 0.261, M is the mass accretion rate onto the surface of the neutron star (equals in 
the young neutron star case), Ct = 1.4 x 10 48 ergs -1 is twice the specific heat at a temperature of 10 9 K, and 
L v9 = 7.4 x 10 39 ergs" 1 is the neutrino cooling rate at T = 10 9 i"C. The terms involving M describe the torque 
and heating by accretion (Levin 1999), and e a is the fraction of accreted energy that goes into heating the star, 
we take it to be unity. The mode energy E = 0.5a 2 tt 2 MR 2 J with J ~ 0.016. The initial LMXB temperature 
is determined self-consistently using equation ( |C3| ) setting a = and T = 0. It varies by a factor of a few over 
three orders of magnitude in accretion rates. The various time-scales in the above equations are listed below. 

Following Lindblom et al. (1999), we write the energy input rate due to gravitational wave back-reaction to 
be 

_L = = _L j,6 (C4) 

with T gr = —18.7s and i^Uz = zAj/lkHz. The negative sign indicates driving. Similarly, the bulk viscosity and 
molecular shear viscosity rates are 

^ = ^vlmTt and i- = ±T 8 - 2 , (C5) 

' bv ' bv ' sv ' sv 

with fbv = 3.5 x 10 17 s and f sv = 2.5 x 10 6 s. We will continue to use these estimates in the presence of a crust 
even though the r-mode occupies a smaller volume. 

BU and this article have pointed to the importance of the boundary layers at the core-crust interface. We 
model the dissipation rate in this layer as a function that is continuous when the layer changes from laminar 
(viscous) to turbulent, 

— = 0(a eq -a) h6(a-a eq ) , (C6) 

T~U T vbl T tbl 
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with 0(x) = 1 when x > and otherwise. The viscous and the turbulent damping rates are, respectively, 

= —V 2 4hz T 8 1 and = — ^CnavkH*, (C7) 

Tvbl Tvbl ' Ttbl T t bl 

where the fiducial time-scales f v bi = 65 s and ftbi = 3.3 x lCP 4 s. We equate r v bi with rtbi to obtain the swith- 

over amplitude a eq = 1.0 x 10~ 3 r)^ 1 T 8 _1 . In our calculations, we evaluate Cd using the Reynolds number 
at the equator. 

The total rate of viscous damping for the r-mode as well as the rate of heating due to the r-mode (ignoring 
bulk viscosity heating) are given by 



